##### Replication for Weidong, Zhang. "The importance of money and connections: 
# Explaining international status from UNGA draft sponsorship networks" (2)

#### Feb, 2025


### last updated  2/12, 2025


setwd("~/Desktop/submission/FDI Networks/International Interaction/Third Revision-1:31:25/Replication 2")
# Note: change it to your own working directory
# Estimated simulation time: 3 hours


# The analysis is conducted by using R version 4.1.1 (2021-08-10) -- "Kick Things"


library("writexl")
library(plyr)
library(dplyr)
library(tidyverse)
library(haven)
library(countrycode)
library(igraph)
library(ergm)
library(btergm) 
#library(janitor)
library(readxl)
library(tibble)
library(ggplot2)
library(texreg)




# Prepare for network in each year====
# Network 2018====
Edge2018 <- read.csv("Edge2018.csv", header=T, row.names=1, stringsAsFactors = FALSE)

# For node attributes (GDP, democracy, military expenditure) of 2018
load("vcov3.RData")

attri2018 <- subset(vcov3, year==2018,
                    select=state: log_military)


ally2018 <- read.csv("ally2018.csv", header=T, row.names=1)
ally2018[,1:196] <- lapply(ally2018[,1:196], as.numeric)

FDI2018 <- read.csv("FDI2018.csv", header=T, row.names=1)
FDI2018[,1:196] <- lapply(FDI2018[,1:196], as.numeric)

trade2018 <- read.csv("trade2018.csv", header=T, row.names=1)
trade2018[,1:196] <- lapply(trade2018[,1:196], as.numeric)


dis2018 <- read.csv("dis2018.csv", header=T, row.names=1)
dis2018[,1:196] <- lapply(dis2018[,1:196], as.numeric)

#set attributes

net2018 <- network(as.matrix(Edge2018), directed = T)

set.vertex.attribute(net2018,names(attri2018),attri2018)


set.network.attribute(net2018,"ally2018",as.matrix(ally2018))

set.network.attribute(net2018,"FDI2018",as.matrix(FDI2018))

set.network.attribute(net2018,"trade2018",as.matrix(trade2018))

set.network.attribute(net2018,"dis2018",as.matrix(dis2018))

# Network 2017====
Edge2017 <- read.csv("Edge2017.csv", header=T, row.names=1)

ally2017 <- read.csv("ally2017.csv", header=T, row.names=1)
ally2017[,1:196] <- lapply(ally2017[,1:196], as.numeric)

FDI2017 <- read.csv("FDI2017.csv", header=T, row.names=1)
FDI2017[,1:196] <- lapply(FDI2017[,1:196], as.numeric)

trade2017 <- read.csv("trade2017.csv", header=T, row.names=1)
trade2017[,1:196] <- lapply(trade2017[,1:196], as.numeric)


dis2017 <- read.csv("dis2017.csv", header=T, row.names=1)
dis2017[,1:196] <- lapply(dis2017[,1:196], as.numeric)


# For node attributes (GDP, democracy, military expenditure) of 2017
load("vcov3.RData")
attri2017 <- subset(vcov3, year==2017,
                    select=state: log_military)



net2017 <- network(as.matrix(Edge2017), directed = T)

set.vertex.attribute(net2017,names(attri2017),attri2017)

set.network.attribute(net2017,"ally2017",as.matrix(ally2017))

set.network.attribute(net2017,"FDI2017",as.matrix(FDI2017))

set.network.attribute(net2017,"trade2017",as.matrix(trade2017))

set.network.attribute(net2017,"dis2017",as.matrix(dis2017))


# Network 2016====
Edge2016 <- read.csv("Edge2016.csv", header=T, row.names=1)
ally2016 <- read.csv("ally2016.csv", header=T, row.names=1)
ally2016[,1:196] <- lapply(ally2016[,1:196], as.numeric)

FDI2016 <- read.csv("FDI2016.csv", header=T, row.names=1)
FDI2016[,1:196] <- lapply(FDI2016[,1:196], as.numeric)

trade2016 <- read.csv("trade2016.csv", header=T, row.names=1)
trade2016[,1:196] <- lapply(trade2016[,1:196], as.numeric)

dis2016 <- read.csv("dis2016.csv", header=T, row.names=1)
dis2016[,1:196] <- lapply(dis2016[,1:196], as.numeric)


# For node attributes (GDP, democracy, military expenditure) of 2016
load("vcov3.RData")
attri2016 <- subset(vcov3, year==2016,
                    select=state: log_military)

net2016 <- network(as.matrix(Edge2016), directed = T)

set.vertex.attribute(net2016,names(attri2016),attri2016)

set.network.attribute(net2016,"ally2016",as.matrix(ally2016))

set.network.attribute(net2016,"FDI2016",as.matrix(FDI2016))

set.network.attribute(net2016,"trade2016",as.matrix(trade2016))

set.network.attribute(net2016,"dis2016",as.matrix(dis2016))

# Network 2015====
Edge2015 <- read.csv("Edge2015.csv", header=T, row.names=1)
ally2015 <- read.csv("ally2015.csv", header=T, row.names=1)
ally2015[,1:196] <- lapply(ally2015[,1:196], as.numeric)

FDI2015 <- read.csv("FDI2015.csv", header=T, row.names=1)
FDI2015[,1:196] <- lapply(FDI2015[,1:196], as.numeric)

trade2015 <- read.csv("trade2015.csv", header=T, row.names=1)
trade2015[,1:196] <- lapply(trade2015[,1:196], as.numeric)

dis2015 <- read.csv("dis2015.csv", header=T, row.names=1)
dis2015[,1:196] <- lapply(dis2015[,1:196], as.numeric)


# For node attributes (GDP, democracy, military expenditure) of 2015
load("vcov3.RData")
attri2015 <- subset(vcov3, year==2015,
                    select=state: log_military)

net2015 <- network(as.matrix(Edge2015), directed = T)

set.vertex.attribute(net2015,names(attri2015),attri2015)

set.network.attribute(net2015,"ally2015",as.matrix(ally2015))

set.network.attribute(net2015,"FDI2015",as.matrix(FDI2015))

set.network.attribute(net2015,"trade2015",as.matrix(trade2015))

set.network.attribute(net2015,"dis2015",as.matrix(dis2015))


# Network 2014====
Edge2014 <- read.csv("Edge2014.csv", header=T, row.names=1)
ally2014 <- read.csv("ally2014.csv", header=T, row.names=1)
ally2014[,1:196] <- lapply(ally2014[,1:196], as.numeric)

FDI2014 <- read.csv("FDI2014.csv", header=T, row.names=1)
FDI2014[,1:196] <- lapply(FDI2014[,1:196], as.numeric)

trade2014 <- read.csv("trade2014.csv", header=T, row.names=1)
trade2014[,1:196] <- lapply(trade2014[,1:196], as.numeric)

dis2014 <- read.csv("dis2014.csv", header=T, row.names=1)
dis2014[,1:196] <- lapply(dis2014[,1:196], as.numeric)


# For node attributes (GDP, democracy, military expenditure) of 2014
load("vcov3.RData")
attri2014 <- subset(vcov3, year==2014,
                    select=state: log_military)

net2014 <- network(as.matrix(Edge2014), directed = T)

set.vertex.attribute(net2014,names(attri2014),attri2014)

set.network.attribute(net2014,"ally2014",as.matrix(ally2014))

set.network.attribute(net2014,"FDI2014",as.matrix(FDI2014))

set.network.attribute(net2014,"trade2014",as.matrix(trade2014))

set.network.attribute(net2014,"dis2014",as.matrix(dis2014))

# Network 2013====
Edge2013 <- read.csv("Edge2013.csv", header=T, row.names=1)
ally2013 <- read.csv("ally2013.csv", header=T, row.names=1)
ally2013[,1:196] <- lapply(ally2013[,1:196], as.numeric)

FDI2013 <- read.csv("FDI2013.csv", header=T, row.names=1)
FDI2013[,1:196] <- lapply(FDI2013[,1:196], as.numeric)
trade2013 <- read.csv("trade2013.csv", header=T, row.names=1)
trade2013[,1:196] <- lapply(trade2013[,1:196], as.numeric)

dis2013 <- read.csv("dis2013.csv", header=T, row.names=1)
dis2013[,1:196] <- lapply(dis2013[,1:196], as.numeric)


# For node attributes (GDP, democracy, military expenditure) of 2013
load("vcov3.RData")
attri2013 <- subset(vcov3, year==2013,
                    select=state: log_military)

net2013 <- network(as.matrix(Edge2013), directed = T)

set.vertex.attribute(net2013,names(attri2013),attri2013)

set.network.attribute(net2013,"ally2013",as.matrix(ally2013))

set.network.attribute(net2013,"FDI2013",as.matrix(FDI2013))
set.network.attribute(net2013,"trade2013",as.matrix(trade2013))

set.network.attribute(net2013,"dis2013",as.matrix(dis2013))


# Network 2012====
Edge2012 <- read.csv("Edge2012.csv", header=T, row.names=1)
ally2012 <- read.csv("ally2012.csv", header=T, row.names=1)
ally2012[,1:196] <- lapply(ally2012[,1:196], as.numeric)

FDI2012 <- read.csv("FDI2012.csv", header=T, row.names=1)
FDI2012[,1:196] <- lapply(FDI2012[,1:196], as.numeric)
trade2012 <- read.csv("trade2012.csv", header=T, row.names=1)
trade2012[,1:196] <- lapply(trade2012[,1:196], as.numeric)

dis2012 <- read.csv("dis2012.csv", header=T, row.names=1)
dis2012[,1:196] <- lapply(dis2012[,1:196], as.numeric)



# For node attributes (GDP, democracy, military expenditure) of 2012
load("vcov3.RData")
attri2012 <- subset(vcov3, year==2012,
                    select=state: log_military)

net2012 <- network(as.matrix(Edge2012), directed = T)

set.vertex.attribute(net2012,names(attri2012),attri2012)

set.network.attribute(net2012,"ally2012",as.matrix(ally2012))

set.network.attribute(net2012,"FDI2012",as.matrix(FDI2012))

set.network.attribute(net2012,"trade2012",as.matrix(trade2012))

set.network.attribute(net2012,"dis2012",as.matrix(dis2012))


# Network 2011====
Edge2011 <- read.csv("Edge2011.csv", header=T, row.names=1)
ally2011 <- read.csv("ally2011.csv", header=T, row.names=1)
ally2011[,1:196] <- lapply(ally2011[,1:196], as.numeric)

FDI2011 <- read.csv("FDI2011.csv", header=T, row.names=1)
FDI2011[,1:196] <- lapply(FDI2011[,1:196], as.numeric)
trade2011 <- read.csv("trade2011.csv", header=T, row.names=1)
trade2011[,1:196] <- lapply(trade2011[,1:196], as.numeric)

dis2011 <- read.csv("dis2011.csv", header=T, row.names=1)
dis2011[,1:196] <- lapply(dis2011[,1:196], as.numeric)


# For node attributes (GDP, democracy, military expenditure) of 2011
load("vcov3.RData")
attri2011 <- subset(vcov3, year==2011,
                    select=state: log_military)

net2011 <- network(as.matrix(Edge2011), directed = T)

set.vertex.attribute(net2011,names(attri2011),attri2011)

set.network.attribute(net2011,"ally2011",as.matrix(ally2011))

set.network.attribute(net2011,"FDI2011",as.matrix(FDI2011))
set.network.attribute(net2011,"trade2011",as.matrix(trade2011))

set.network.attribute(net2011,"dis2011",as.matrix(dis2011))


# Network 2010====
Edge2010 <- read.csv("Edge2010.csv", header=T, row.names=1)
ally2010 <- read.csv("ally2010.csv", header=T, row.names=1)
ally2010[,1:196] <- lapply(ally2010[,1:196], as.numeric)

FDI2010 <- read.csv("FDI2010.csv", header=T, row.names=1)
FDI2010[,1:196] <- lapply(FDI2010[,1:196], as.numeric)
trade2010 <- read.csv("trade2010.csv", header=T, row.names=1)
trade2010[,1:196] <- lapply(trade2010[,1:196], as.numeric)

dis2010 <- read.csv("dis2010.csv", header=T, row.names=1)
dis2010[,1:196] <- lapply(dis2010[,1:196], as.numeric)


# For node attributes (GDP, democracy, military expenditure) of 2010
load("vcov3.RData")
attri2010 <- subset(vcov3, year==2010,
                    select=state: log_military)

net2010 <- network(as.matrix(Edge2010), directed = T)

set.vertex.attribute(net2010,names(attri2010),attri2010)

set.network.attribute(net2010,"ally2010",as.matrix(ally2010))

set.network.attribute(net2010,"FDI2010",as.matrix(FDI2010))
set.network.attribute(net2010,"trade2010",as.matrix(trade2010))

set.network.attribute(net2010,"dis2010",as.matrix(dis2010))


# Network 2009====
Edge2009 <- read.csv("Edge2009.csv", header=T, row.names=1)
ally2009 <- read.csv("ally2009.csv", header=T, row.names=1)
ally2009[,1:196] <- lapply(ally2009[,1:196], as.numeric)

FDI2009 <- read.csv("FDI2009.csv", header=T, row.names=1)
FDI2009[,1:196] <- lapply(FDI2009[,1:196], as.numeric)
trade2009 <- read.csv("trade2009.csv", header=T, row.names=1)
trade2009[,1:196] <- lapply(trade2009[,1:196], as.numeric)

dis2009 <- read.csv("dis2009.csv", header=T, row.names=1)
dis2009[,1:196] <- lapply(dis2009[,1:196], as.numeric)


# For node attributes (GDP, democracy, military expenditure) of 2009
load("vcov3.RData")
attri2009 <- subset(vcov3, year==2009,
                    select=state: log_military)

net2009 <- network(as.matrix(Edge2009), directed = T)

set.vertex.attribute(net2009,names(attri2009),attri2009)

set.network.attribute(net2009,"ally2009",as.matrix(ally2009))

set.network.attribute(net2009,"FDI2009",as.matrix(FDI2009))
set.network.attribute(net2009,"trade2009",as.matrix(trade2009))

set.network.attribute(net2009,"dis2009",as.matrix(dis2009))



##### Now Using BTERGM to do the same thing
networks <- list()# create network object
networks[[1]] <- net2010
networks[[2]] <- net2011
networks[[3]] <- net2012
networks[[4]] <- net2013
networks[[5]] <- net2014
networks[[6]] <- net2015
networks[[7]] <- net2016
networks[[8]] <- net2017
networks[[9]] <- net2018

# Set network attributes
network.vertex.names(x=net2009) <- attri2009$state
network.vertex.names(x=net2010) <- attri2010$state
network.vertex.names(x=net2011) <- attri2011$state
network.vertex.names(x=net2012) <- attri2012$state
network.vertex.names(x=net2013) <- attri2013$state
network.vertex.names(x=net2014) <- attri2014$state
network.vertex.names(x=net2015) <- attri2015$state
network.vertex.names(x=net2016) <- attri2016$state
network.vertex.names(x=net2017) <- attri2017$state
network.vertex.names(x=net2018) <- attri2018$state


LagNet<-list()
LagNet[[1]] <- net2009      # add network to the list
LagNet[[2]] <- net2010
LagNet[[3]] <- net2011
LagNet[[4]] <- net2012
LagNet[[5]] <- net2013
LagNet[[6]] <- net2014
LagNet[[7]] <- net2015
LagNet[[8]] <- net2016
LagNet[[9]] <- net2017


FDI2009<-network(as.matrix(FDI2009),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight", directed = T)

FDI2010<-network(as.matrix(FDI2010),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight", directed = T)


FDI2011<-network(as.matrix(FDI2011),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight", directed = T)

FDI2012<-network(as.matrix(FDI2012),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight", directed = T)

FDI2013<-network(as.matrix(FDI2013),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight", directed = T)

FDI2014<-network(as.matrix(FDI2014),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight", directed = T)

FDI2015<-network(as.matrix(FDI2015),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight", directed = T)

FDI2016<-network(as.matrix(FDI2016),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight", directed = T)

FDI2017<-network(as.matrix(FDI2017),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight", directed = T)

FDI2018<-network(as.matrix(FDI2018),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight", directed = T)


# Set network attributes
network.vertex.names(x=FDI2009) <- attri2009$state
network.vertex.names(x=FDI2010) <- attri2010$state
network.vertex.names(x=FDI2011) <- attri2011$state
network.vertex.names(x=FDI2012) <- attri2012$state
network.vertex.names(x=FDI2013) <- attri2013$state
network.vertex.names(x=FDI2014) <- attri2014$state
network.vertex.names(x=FDI2015) <- attri2015$state
network.vertex.names(x=FDI2016) <- attri2016$state
network.vertex.names(x=FDI2017) <- attri2017$state
network.vertex.names(x=FDI2018) <- attri2018$state


FDINet<-list()

FDINet[[1]] <- FDI2009
FDINet[[2]] <- FDI2010
FDINet[[3]] <- FDI2011
FDINet[[4]] <- FDI2012
FDINet[[5]] <- FDI2013
FDINet[[6]] <- FDI2014
FDINet[[7]] <- FDI2015
FDINet[[8]] <- FDI2016
FDINet[[9]] <- FDI2017

# For ally
ally2009<-network(as.matrix(ally2009), directed = F)
ally2010<-network(as.matrix(ally2010), directed = F)
ally2011<-network(as.matrix(ally2011), directed = F)
ally2012<-network(as.matrix(ally2012), directed = F)
ally2013<-network(as.matrix(ally2013), directed = F)
ally2014<-network(as.matrix(ally2014), directed = F)
ally2015<-network(as.matrix(ally2015), directed = F)
ally2016<-network(as.matrix(ally2016), directed = F)
ally2017<-network(as.matrix(ally2017), directed = F)
ally2018<-network(as.matrix(ally2018), directed = F)

# Set network attributes
network.vertex.names(x=ally2009) <- attri2009$state
network.vertex.names(x=ally2010) <- attri2010$state
network.vertex.names(x=ally2011) <- attri2011$state
network.vertex.names(x=ally2012) <- attri2012$state
network.vertex.names(x=ally2013) <- attri2013$state
network.vertex.names(x=ally2014) <- attri2014$state
network.vertex.names(x=ally2015) <- attri2015$state
network.vertex.names(x=ally2016) <- attri2016$state
network.vertex.names(x=ally2017) <- attri2017$state
network.vertex.names(x=ally2018) <- attri2018$state




allyNet<-list()
allyNet[[1]] <-ally2010
allyNet[[2]] <-ally2011
allyNet[[3]] <-ally2012
allyNet[[4]] <-ally2013
allyNet[[5]] <-ally2014
allyNet[[6]] <-ally2015
allyNet[[7]] <-ally2016
allyNet[[8]] <-ally2017
allyNet[[9]] <-ally2018


# For trade
trade2009<-network(as.matrix(trade2009),matrix.type="a", 
                   ignore.eval=FALSE,names.eval="weight_trade", directed = T)

trade2010<-network(as.matrix(trade2010),matrix.type="a", 
                   ignore.eval=FALSE,names.eval="weight_trade", directed = T)


trade2011<-network(as.matrix(trade2011),matrix.type="a", 
                   ignore.eval=FALSE,names.eval="weight_trade", directed = T)

trade2012<-network(as.matrix(trade2012),matrix.type="a", 
                   ignore.eval=FALSE,names.eval="weight_trade", directed = T)

trade2013<-network(as.matrix(trade2013),matrix.type="a", 
                   ignore.eval=FALSE,names.eval="weight_trade", directed = T)

trade2014<-network(as.matrix(trade2014),matrix.type="a", 
                   ignore.eval=FALSE,names.eval="weight_trade", directed = T)

trade2015<-network(as.matrix(trade2015),matrix.type="a", 
                   ignore.eval=FALSE,names.eval="weight_trade", directed = T)

trade2016<-network(as.matrix(trade2016),matrix.type="a", 
                   ignore.eval=FALSE,names.eval="weight_trade", directed = T)

trade2017<-network(as.matrix(trade2017),matrix.type="a", 
                   ignore.eval=FALSE,names.eval="weight_trade", directed = T)

trade2018<-network(as.matrix(trade2018),matrix.type="a", 
                   ignore.eval=FALSE,names.eval="weight_trade", directed = T)



# Set network attributes
network.vertex.names(x=trade2009) <- attri2009$state
network.vertex.names(x=trade2010) <- attri2010$state
network.vertex.names(x=trade2011) <- attri2011$state
network.vertex.names(x=trade2012) <- attri2012$state
network.vertex.names(x=trade2013) <- attri2013$state
network.vertex.names(x=trade2014) <- attri2014$state
network.vertex.names(x=trade2015) <- attri2015$state
network.vertex.names(x=trade2016) <- attri2016$state
network.vertex.names(x=trade2017) <- attri2017$state
network.vertex.names(x=trade2018) <- attri2018$state

# 
tradeNet<-list()

tradeNet[[1]] <- trade2009
tradeNet[[2]] <- trade2010
tradeNet[[3]] <- trade2011
tradeNet[[4]] <- trade2012
tradeNet[[5]] <- trade2013
tradeNet[[6]] <- trade2014
tradeNet[[7]] <- trade2015
tradeNet[[8]] <- trade2016
tradeNet[[9]] <- trade2017

# For distance
dis2009<-network(as.matrix(dis2009),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight_distance", directed = T)

dis2010<-network(as.matrix(dis2010),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight_distance", directed = T)


dis2011<-network(as.matrix(dis2011),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight_distance", directed = T)

dis2012<-network(as.matrix(dis2012),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight_distance", directed = T)

dis2013<-network(as.matrix(dis2013),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight_distance", directed = T)

dis2014<-network(as.matrix(dis2014),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight_distance", directed = T)

dis2015<-network(as.matrix(dis2015),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight_distance", directed = T)

dis2016<-network(as.matrix(dis2016),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight_distance", directed = T)

dis2017<-network(as.matrix(dis2017),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight_distance", directed = T)

dis2018<-network(as.matrix(dis2018),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight_distance", directed = T)

# Set network attributes
network.vertex.names(x=dis2009) <- attri2009$state
network.vertex.names(x=dis2010) <- attri2010$state
network.vertex.names(x=dis2011) <- attri2011$state
network.vertex.names(x=dis2012) <- attri2012$state
network.vertex.names(x=dis2013) <- attri2013$state
network.vertex.names(x=dis2014) <- attri2014$state
network.vertex.names(x=dis2015) <- attri2015$state
network.vertex.names(x=dis2016) <- attri2016$state
network.vertex.names(x=dis2017) <- attri2017$state
network.vertex.names(x=dis2018) <- attri2018$state

disNet<-list()
disNet[[1]] <-dis2010
disNet[[2]] <-dis2011
disNet[[3]] <-dis2012
disNet[[4]] <-dis2013
disNet[[5]] <-dis2014
disNet[[6]] <-dis2015
disNet[[7]] <-dis2016
disNet[[8]] <-dis2017
disNet[[9]] <-dis2018


# set overall network attributes

set.vertex.attribute(networks[[1]],"attri", attri2010)

set.vertex.attribute(networks[[2]],"attri", attri2011)

set.vertex.attribute(networks[[3]],"attri", attri2012)
set.vertex.attribute(networks[[4]],"attri", attri2013)
set.vertex.attribute(networks[[5]],"attri", attri2014)
set.vertex.attribute(networks[[6]],"attri", attri2015)
set.vertex.attribute(networks[[7]],"attri", attri2016)
set.vertex.attribute(networks[[8]],"attri", attri2017)
set.vertex.attribute(networks[[9]],"attri", attri2018)



# Model 2 for Table 2 ====
set.seed(1234)
#
fit2 <- btergm(networks ~
                 edgecov(FDINet, attrname="weight")
               +edgecov(LagNet)
               +edgecov(tradeNet,attrname="weight_trade")
               +edgecov(disNet, attrname="weight_distance")
               
               +edgecov(allyNet)
               
               + edges+ istar(2) + ostar(2) + mutual + ttriple + ctriple
               
               +nodeicov("log_military")
               +nodeicov("mean_Democracy")
               +nodeocov("log_GDP")
               
               +absdiff("mean_Democracy")
               +absdiff("log_GDP")
               +absdiff("log_military")
               +isolates, 
               R = 1000)
summary(fit2)
# Note:"Boot mean", not "Estimate", is used in Table 2 of the paper



################################
# Get the FDI graph for model 2 in Table 2: Lower left plot of Figure 3====

alldyads2 <- edgeprob(fit2)

model <- glm(probability ~ edgecov.weight, family = "binomial", data = alldyads2)

# Create a data frame for predictions
new_data <- alldyads2 %>%
  mutate(pred = predict(model, newdata = ., type = "response"),
         se = predict(model, newdata = ., type = "link", se.fit = TRUE)$se.fit) %>%
  mutate(lower = pred - 1.96 * se,  # Lower bound of 95% CI
         upper = pred + 1.96 * se)  # Upper bound of 95% CI

# Plot with ggplot2
ggplot(new_data, aes(x = edgecov.weight, y = pred)) +
  geom_line(color = "blue") +                                # Main line
  geom_line(aes(y = lower), linetype = "dashed", color = "black") +  # Lower bound
  geom_line(aes(y = upper), linetype = "dashed", color = "black") +  # Upper bound
  labs(y = "Probability of Cosponsorship", x = "FDI (log)",caption = 
         "The effect of FDI of Model 2") +
  theme_minimal()+
  theme(
    axis.title.x = element_text(size = 24),  # X-axis label size
    axis.title.y = element_text(size = 24),  # Y-axis label size
    axis.text.x = element_text(size = 24),   # X-axis numbers (ticks) size
    axis.text.y = element_text(size = 24),    # Y-axis numbers (ticks) size
    plot.caption = element_text(size = 24, face = "bold", hjust = 0.5)  # Bigger, centered caption
  )


# Get the effect of trade: Lower right plot of Figure 3 ====

model <- glm(probability ~ edgecov.weight_trade, family = "binomial", data = alldyads2)

# Create a data frame for predictions
new_data <- alldyads2 %>%
  mutate(pred = predict(model, newdata = ., type = "response"),
         se = predict(model, newdata = ., type = "link", se.fit = TRUE)$se.fit) %>%
  mutate(lower = pred - 1.96 * se,  # Lower bound of 95% CI
         upper = pred + 1.96 * se)  # Upper bound of 95% CI

# Plot with ggplot2
ggplot(new_data, aes(x = edgecov.weight_trade, y = pred)) +
  geom_line(color = "blue") +                                # Main line
  geom_line(aes(y = lower), linetype = "dashed", color = "black") +  # Lower bound
  geom_line(aes(y = upper), linetype = "dashed", color = "black") +  # Upper bound
  labs(y = "Probability of Cosponsorship", x = "Trade (log)",caption = 
         "The effect of trade of Model 2") +
  theme_minimal()+
  theme(
    axis.title.x = element_text(size = 24),  # X-axis label size
    axis.title.y = element_text(size = 24),  # Y-axis label size
    axis.text.x = element_text(size = 24),   # X-axis numbers (ticks) size
    axis.text.y = element_text(size = 24),    # Y-axis numbers (ticks) size
    plot.caption = element_text(size = 24, face = "bold", hjust = 0.5)  # Bigger, centered caption
  )

# goodness of fit: for Figure A2 ====

gof2<-gof(fit2, control=control.gof.ergm(nsim=1000))
plot(gof2)


# End of script



























